T    = 1.0d2
y    = 0.0d0
nH   = 3.7771377550788764E-3
G0   = 1.5891954714728652E-11
tau0 = 555.0763465450948

sigma0 = 6.3d-18

N=1000

RC = hiirca_hui(T)
CI = hici_hui(T)

C =  G0 / (RC * nH)

xHII_0 = (-C + sqrt(C*C + 4*C) ) * 0.5d0


; want x (ionized fraction) to go from xHII_0 to small

xHII_min = 1.0d-15
xHII_ana = dblarr(N)

xHII_ana = xHII_0 - dindgen(N)/(N-1) * (xHII_0 - xHII_min) 

lhs = 2 * alog( (1.0d0-xHII_ana) / xHII_ana  * xHII_0 / (1.0d0-xHII_0) ) - $
      ( 1.0d0 / (1.0d0-xHII_ana) - 1 / (1.0d0-xHII_0) ) - tau0

lhs0 = 2 * alog( (1.0d0-xHII_ana) / xHII_ana  * xHII_0 / (1.0d0-xHII_0) ) - $
      ( 1.0d0 / (1.0d0-xHII_ana) - 1 / (1.0d0-xHII_0) )


altay_set_x
altay_plot, alog10(xHII_ana), lhs

altay_oplot, [-100,100], [0,0], color=100
altay_oplot, [-10,-10], [-1000,1000], color=100










end
